home *** CD-ROM | disk | FTP | other *** search
/ Disc to the Future 2 / Disc to the Future Part II Programmer's Reference (Wayzata Technology)(6013)(1992).bin / UNIX / GNU / GNU-PLOT / STANDARD.C < prev    next >
C/C++ Source or Header  |  1992-01-06  |  16KB  |  803 lines
  1. /* GNUPLOT - standard.c */
  2. /*
  3.  * Copyright (C) 1986, 1987, 1990, 1991   Thomas Williams, Colin Kelley
  4.  *
  5.  * Permission to use, copy, and distribute this software and its
  6.  * documentation for any purpose with or without fee is hereby granted, 
  7.  * provided that the above copyright notice appear in all copies and 
  8.  * that both that copyright notice and this permission notice appear 
  9.  * in supporting documentation.
  10.  *
  11.  * Permission to modify the software is granted, but not the right to
  12.  * distribute the modified code.  Modifications are to be distributed 
  13.  * as patches to released version.
  14.  *  
  15.  * This software is provided "as is" without express or implied warranty.
  16.  * 
  17.  *
  18.  * AUTHORS
  19.  * 
  20.  *   Original Software:
  21.  *     Thomas Williams,  Colin Kelley.
  22.  * 
  23.  *   Gnuplot 2.0 additions:
  24.  *       Russell Lang, Dave Kotz, John Campbell.
  25.  *
  26.  *   Gnuplot 3.0 additions:
  27.  *       Gershon Elber and many others.
  28.  * 
  29.  * Send your comments or suggestions to 
  30.  *  pixar!info-gnuplot@sun.com.
  31.  * This is a mailing list; to join it send a note to 
  32.  *  pixar!info-gnuplot-request@sun.com.  
  33.  * Send bug reports to
  34.  *  pixar!bug-gnuplot@sun.com.
  35.  */
  36.  
  37. #include <math.h>
  38. #include <stdio.h>
  39. #include "plot.h"
  40.  
  41. #ifdef vms
  42. #include <errno.h>
  43. #else
  44. extern int errno;
  45. #endif /* vms */
  46.  
  47.  
  48. extern struct value stack[STACK_DEPTH];
  49. extern int s_p;
  50. extern double zero;
  51.  
  52. struct value *pop(), *complex(), *integer();
  53.  
  54. double magnitude(), angle(), real(), imag();
  55.  
  56. /* The bessel function approximations here are from
  57.  * "Computer Approximations"
  58.  * by Hart, Cheney et al.
  59.  * John Wiley & Sons, 1968
  60.  */
  61.  
  62. /* There appears to be a mistake in Hart, Cheney et al. on page 149.
  63.  * Where it list Qn(x)/x ~ P(z*z)/Q(z*z), z = 8/x, it should read
  64.  *               Qn(x)/z ~ P(z*z)/Q(z*z), z = 8/x
  65.  * In the functions below, Qn(x) is implementated using the later
  66.  * equation.
  67.  * These bessel functions are accurate to about 1e-13
  68.  */
  69.  
  70. #define PI_ON_FOUR       0.78539816339744830961566084581987572
  71. #define PI_ON_TWO        1.57079632679489661923131269163975144
  72. #define THREE_PI_ON_FOUR 2.35619449019234492884698253745962716
  73. #define TWO_ON_PI        0.63661977236758134307553505349005744
  74.  
  75. static double dzero = 0.0;
  76.  
  77. /* jzero for x in [0,8]
  78.  * Index 5849, 19.22 digits precision
  79.  */
  80. static double pjzero[] = {
  81.      0.4933787251794133561816813446e+21,
  82.     -0.11791576291076105360384408e+21,
  83.      0.6382059341072356562289432465e+19,
  84.     -0.1367620353088171386865416609e+18,
  85.      0.1434354939140346111664316553e+16,
  86.     -0.8085222034853793871199468171e+13,
  87.      0.2507158285536881945555156435e+11,
  88.     -0.4050412371833132706360663322e+8,
  89.      0.2685786856980014981415848441e+5
  90. };
  91.  
  92. static double qjzero[] = {
  93.     0.4933787251794133562113278438e+21,
  94.     0.5428918384092285160200195092e+19,
  95.     0.3024635616709462698627330784e+17,
  96.     0.1127756739679798507056031594e+15,
  97.     0.3123043114941213172572469442e+12,
  98.     0.669998767298223967181402866e+9,
  99.     0.1114636098462985378182402543e+7,
  100.     0.1363063652328970604442810507e+4,
  101.     0.1e+1
  102. };
  103.  
  104. /* pzero for x in [8,inf]
  105.  * Index 6548, 18.16 digits precision
  106.  */
  107. static double ppzero[] = {
  108.     0.2277909019730468430227002627e+5,
  109.     0.4134538663958076579678016384e+5,
  110.     0.2117052338086494432193395727e+5,
  111.     0.348064864432492703474453111e+4,
  112.     0.15376201909008354295771715e+3,
  113.     0.889615484242104552360748e+0
  114. };
  115.  
  116. static double qpzero[] = {
  117.     0.2277909019730468431768423768e+5,
  118.     0.4137041249551041663989198384e+5,
  119.     0.2121535056188011573042256764e+5,
  120.     0.350287351382356082073561423e+4,
  121.     0.15711159858080893649068482e+3,
  122.     0.1e+1
  123. };
  124.  
  125. /* qzero for x in [8,inf]
  126.  * Index 6948, 18.33 digits precision
  127.  */
  128. static double pqzero[] = {
  129.     -0.8922660020080009409846916e+2,
  130.     -0.18591953644342993800252169e+3,
  131.     -0.11183429920482737611262123e+3,
  132.     -0.2230026166621419847169915e+2,
  133.     -0.124410267458356384591379e+1,
  134.     -0.8803330304868075181663e-2,
  135. };
  136.  
  137. static double qqzero[] = {
  138.     0.571050241285120619052476459e+4,
  139.     0.1195113154343461364695265329e+5,
  140.     0.726427801692110188369134506e+4,
  141.     0.148872312322837565816134698e+4,
  142.     0.9059376959499312585881878e+2,
  143.     0.1e+1
  144. };
  145.  
  146.  
  147. /* yzero for x in [0,8]
  148.  * Index 6245, 18.78 digits precision
  149.  */
  150. static double pyzero[] = {
  151.     -0.2750286678629109583701933175e+20,
  152.      0.6587473275719554925999402049e+20,
  153.     -0.5247065581112764941297350814e+19,
  154.      0.1375624316399344078571335453e+18,
  155.     -0.1648605817185729473122082537e+16,
  156.      0.1025520859686394284509167421e+14,
  157.     -0.3436371222979040378171030138e+11,
  158.      0.5915213465686889654273830069e+8,
  159.     -0.4137035497933148554125235152e+5
  160. };
  161.  
  162. static double qyzero[] = {
  163.     0.3726458838986165881989980739e+21,
  164.     0.4192417043410839973904769661e+19,
  165.     0.2392883043499781857439356652e+17,
  166.     0.9162038034075185262489147968e+14,
  167.     0.2613065755041081249568482092e+12,
  168.     0.5795122640700729537380087915e+9,
  169.     0.1001702641288906265666651753e+7,
  170.     0.1282452772478993804176329391e+4,
  171.     0.1e+1
  172. };
  173.  
  174.  
  175. /* jone for x in [0,8]
  176.  * Index 6050, 20.98 digits precision
  177.  */
  178. static double pjone[] = {
  179.      0.581199354001606143928050809e+21,
  180.     -0.6672106568924916298020941484e+20,
  181.      0.2316433580634002297931815435e+19,
  182.     -0.3588817569910106050743641413e+17,
  183.      0.2908795263834775409737601689e+15,
  184.     -0.1322983480332126453125473247e+13,
  185.      0.3413234182301700539091292655e+10,
  186.     -0.4695753530642995859767162166e+7,
  187.      0.270112271089232341485679099e+4
  188. };
  189.  
  190. static double qjone[] = {
  191.     0.11623987080032122878585294e+22,
  192.     0.1185770712190320999837113348e+20,
  193.     0.6092061398917521746105196863e+17,
  194.     0.2081661221307607351240184229e+15,
  195.     0.5243710262167649715406728642e+12,
  196.     0.1013863514358673989967045588e+10,
  197.     0.1501793594998585505921097578e+7,
  198.     0.1606931573481487801970916749e+4,
  199.     0.1e+1
  200. };
  201.  
  202.  
  203. /* pone for x in [8,inf]
  204.  * Index 6749, 18.11 digits precision
  205.  */
  206. static double ppone[] = {
  207.     0.352246649133679798341724373e+5,
  208.     0.62758845247161281269005675e+5,
  209.     0.313539631109159574238669888e+5,
  210.     0.49854832060594338434500455e+4,
  211.     0.2111529182853962382105718e+3,
  212.     0.12571716929145341558495e+1
  213. };
  214.  
  215. static double qpone[] = {
  216.     0.352246649133679798068390431e+5,
  217.     0.626943469593560511888833731e+5,
  218.     0.312404063819041039923015703e+5,
  219.     0.4930396490181088979386097e+4,
  220.     0.2030775189134759322293574e+3,
  221.     0.1e+1
  222. };
  223.  
  224. /* qone for x in [8,inf]
  225.  * Index 7149, 18.28 digits precision
  226.  */
  227. static double pqone[] = {
  228.     0.3511751914303552822533318e+3,
  229.     0.7210391804904475039280863e+3,
  230.     0.4259873011654442389886993e+3,
  231.     0.831898957673850827325226e+2,
  232.     0.45681716295512267064405e+1,
  233.     0.3532840052740123642735e-1
  234. };
  235.  
  236. static double qqone[] = {
  237.     0.74917374171809127714519505e+4,
  238.     0.154141773392650970499848051e+5,
  239.     0.91522317015169922705904727e+4,
  240.     0.18111867005523513506724158e+4,
  241.     0.1038187585462133728776636e+3,
  242.     0.1e+1
  243. };
  244.  
  245.  
  246. /* yone for x in [0,8]
  247.  * Index 6444, 18.24 digits precision
  248.  */
  249. static double pyone[] = {
  250.     -0.2923821961532962543101048748e+20,
  251.      0.7748520682186839645088094202e+19,
  252.     -0.3441048063084114446185461344e+18,
  253.      0.5915160760490070618496315281e+16,
  254.     -0.4863316942567175074828129117e+14,
  255.      0.2049696673745662182619800495e+12,
  256.     -0.4289471968855248801821819588e+9,
  257.      0.3556924009830526056691325215e+6
  258. };
  259.  
  260. static double qyone[] = {
  261.     0.1491311511302920350174081355e+21,
  262.     0.1818662841706134986885065935e+19,
  263.     0.113163938269888452690508283e+17,
  264.     0.4755173588888137713092774006e+14,
  265.     0.1500221699156708987166369115e+12,
  266.     0.3716660798621930285596927703e+9,
  267.     0.726914730719888456980191315e+6,
  268.     0.10726961437789255233221267e+4,
  269.     0.1e+1
  270. };
  271.  
  272.  
  273. f_real()
  274. {
  275. struct value a;
  276.     push( complex(&a,real(pop(&a)), 0.0) );
  277. }
  278.  
  279. f_imag()
  280. {
  281. struct value a;
  282.     push( complex(&a,imag(pop(&a)), 0.0) );
  283. }
  284.  
  285. f_arg()
  286. {
  287. struct value a;
  288.     push( complex(&a,angle(pop(&a)), 0.0) );
  289. }
  290.  
  291. f_conjg()
  292. {
  293. struct value a;
  294.     (void) pop(&a);
  295.     push( complex(&a,real(&a),-imag(&a) ));
  296. }
  297.  
  298. f_sin()
  299. {
  300. struct value a;
  301.     (void) pop(&a);
  302.     push( complex(&a,sin(real(&a))*cosh(imag(&a)), cos(real(&a))*sinh(imag(&a))) );
  303. }
  304.  
  305. f_cos()
  306. {
  307. struct value a;
  308.     (void) pop(&a);
  309.     push( complex(&a,cos(real(&a))*cosh(imag(&a)), -sin(real(&a))*sinh(imag(&a))));
  310. }
  311.  
  312. f_tan()
  313. {
  314. struct value a;
  315. register double den;
  316.     (void) pop(&a);
  317.     if (imag(&a) == 0.0)
  318.         push( complex(&a,tan(real(&a)),0.0) );
  319.     else {
  320.         den = cos(2*real(&a))+cosh(2*imag(&a));
  321.         if (den == 0.0) {
  322.             undefined = TRUE;
  323.             push( &a );
  324.         }
  325.         else
  326.             push( complex(&a,sin(2*real(&a))/den, sinh(2*imag(&a))/den) );
  327.     }
  328. }
  329.  
  330. f_asin()
  331. {
  332. struct value a;
  333. register double alpha, beta, x, y;
  334.     (void) pop(&a);
  335.     x = real(&a); y = imag(&a);
  336.     if (y == 0.0) {
  337.         if (fabs(x) > 1.0) {
  338.             undefined = TRUE;
  339.             push(complex(&a,0.0, 0.0));
  340.         } else
  341.             push( complex(&a,asin(x),0.0) );
  342.     } else {
  343.         beta  = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
  344.         alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
  345.         push( complex(&a,asin(beta), log(alpha + sqrt(alpha*alpha-1))) );
  346.     }
  347. }
  348.  
  349. f_acos()
  350. {
  351. struct value a;
  352. register double alpha, beta, x, y;
  353.     (void) pop(&a);
  354.     x = real(&a); y = imag(&a);
  355.     if (y == 0.0) {
  356.         if (fabs(x) > 1.0) {
  357.             undefined = TRUE;
  358.             push(complex(&a,0.0, 0.0));
  359.         } else
  360.             push( complex(&a,acos(x),0.0) );
  361.     } else {
  362.         alpha = sqrt((x + 1)*(x + 1) + y*y)/2 + sqrt((x - 1)*(x - 1) + y*y)/2;
  363.         beta  = sqrt((x + 1)*(x + 1) + y*y)/2 - sqrt((x - 1)*(x - 1) + y*y)/2;
  364.         push( complex(&a,acos(beta), log(alpha + sqrt(alpha*alpha-1))) );
  365.     }
  366. }
  367.  
  368. f_atan()
  369. {
  370. struct value a;
  371. register double x, y, u, v, w, z;
  372.     (void) pop(&a);
  373.     x = real(&a); y = imag(&a);
  374.     if (y == 0.0)
  375.         push( complex(&a,atan(x), 0.0) );
  376.     else if (x == 0.0 && fabs(y) == 1.0) {
  377.         undefined = TRUE;
  378.         push(complex(&a,0.0, 0.0));
  379.     } else {
  380.             if (x >= 0) {
  381.                 u = x;
  382.             v = y;
  383.         } else {
  384.                 u = -x;
  385.             v = -y;
  386.         }
  387.         
  388.             z = atan(2*u/(1-u*u-v*v));
  389.         w = log((u*u+(v+1)*(v+1))/(u*u+(v-1)*(v-1)))/4;
  390.         if (z < 0)
  391.                 z = z + 2*PI_ON_TWO;
  392.         if (x < 0) {
  393.                 z = -z;
  394.             w = -w;
  395.         }
  396.         push( complex(&a,0.5*z, w) );
  397.     }
  398. }
  399.  
  400. f_sinh()
  401. {
  402. struct value a;
  403.     (void) pop(&a);
  404.     push( complex(&a,sinh(real(&a))*cos(imag(&a)), cosh(real(&a))*sin(imag(&a))) );
  405. }
  406.  
  407. f_cosh()
  408. {
  409. struct value a;
  410.     (void) pop(&a);
  411.     push( complex(&a,cosh(real(&a))*cos(imag(&a)), sinh(real(&a))*sin(imag(&a))) );
  412. }
  413.  
  414. f_tanh()
  415. {
  416. struct value a;
  417. register double den;
  418.     (void) pop(&a);
  419.     den = cosh(2*real(&a)) + cos(2*imag(&a));
  420.     push( complex(&a,sinh(2*real(&a))/den, sin(2*imag(&a))/den) );
  421. }
  422.  
  423. f_int()
  424. {
  425. struct value a;
  426.     push( integer(&a,(int)real(pop(&a))) );
  427. }
  428.  
  429.  
  430. f_abs()
  431. {
  432. struct value a;
  433.     (void) pop(&a);
  434.     switch (a.type) {
  435.         case INT:
  436.             push( integer(&a,abs(a.v.int_val)) );            
  437.             break;
  438.         case CMPLX:
  439.             push( complex(&a,magnitude(&a), 0.0) );
  440.     }
  441. }
  442.  
  443. f_sgn()
  444. {
  445. struct value a;
  446.     (void) pop(&a);
  447.     switch(a.type) {
  448.         case INT:
  449.             push( integer(&a,(a.v.int_val > 0) ? 1 : 
  450.                     (a.v.int_val < 0) ? -1 : 0) );
  451.             break;
  452.         case CMPLX:
  453.             push( integer(&a,(a.v.cmplx_val.real > 0.0) ? 1 : 
  454.                     (a.v.cmplx_val.real < 0.0) ? -1 : 0) );
  455.             break;
  456.     }
  457. }
  458.  
  459.  
  460. f_sqrt()
  461. {
  462. struct value a;
  463. register double mag, ang;
  464.     (void) pop(&a);
  465.     mag = sqrt(magnitude(&a));
  466.     if (imag(&a) == 0.0 && real(&a) < 0.0)
  467.         push( complex(&a,0.0,mag) );
  468.     else
  469.     {
  470.         if ( (ang = angle(&a)) < 0.0)
  471.             ang += 2*Pi;
  472.         ang /= 2;
  473.         push( complex(&a,mag*cos(ang), mag*sin(ang)) );
  474.     }
  475. }
  476.  
  477.  
  478. f_exp()
  479. {
  480. struct value a;
  481. register double mag, ang;
  482.     (void) pop(&a);
  483.     mag = exp(real(&a));
  484.     ang = imag(&a);
  485.     push( complex(&a,mag*cos(ang), mag*sin(ang)) );
  486. }
  487.  
  488.  
  489. f_log10()
  490. {
  491. struct value a;
  492. register double l10;;
  493.     (void) pop(&a);
  494.     l10 = log(10.0);    /***** replace with a constant! ******/
  495.     push( complex(&a,log(magnitude(&a))/l10, angle(&a)/l10) );
  496. }
  497.  
  498.  
  499. f_log()
  500. {
  501. struct value a;
  502.     (void) pop(&a);
  503.     push( complex(&a,log(magnitude(&a)), angle(&a)) );
  504. }
  505.  
  506.  
  507. f_floor()
  508. {
  509. struct value a;
  510.  
  511.     (void) pop(&a);
  512.     switch (a.type) {
  513.         case INT:
  514.             push( integer(&a,(int)floor((double)a.v.int_val)));            
  515.             break;
  516.         case CMPLX:
  517.             push( integer(&a,(int)floor(a.v.cmplx_val.real)));
  518.     }
  519. }
  520.  
  521.  
  522. f_ceil()
  523. {
  524. struct value a;
  525.  
  526.     (void) pop(&a);
  527.     switch (a.type) {
  528.         case INT:
  529.             push( integer(&a,(int)ceil((double)a.v.int_val)));            
  530.             break;
  531.         case CMPLX:
  532.             push( integer(&a,(int)ceil(a.v.cmplx_val.real)));
  533.     }
  534. }
  535.  
  536. #ifdef GAMMA
  537.  
  538. f_gamma()
  539. {
  540. extern int signgam;
  541. register double y;
  542. struct value a;
  543.  
  544.     y = GAMMA(real(pop(&a)));
  545.     if (y > 88.0) {
  546.         undefined = TRUE;
  547.         push( integer(&a,0) );
  548.     }
  549.     else
  550.         push( complex(&a,signgam * exp(y),0.0) );
  551. }
  552.  
  553. #endif /* GAMMA */
  554.  
  555.  
  556. /* bessel function approximations */
  557. double jzero(x)
  558. double x;
  559. {
  560. double p, q, x2;
  561. int n;
  562.  
  563.     x2 = x * x;
  564.     p = pjzero[8];
  565.     q = qjzero[8];
  566.     for (n=7; n>=0; n--) {
  567.         p = p*x2 + pjzero[n];
  568.         q = q*x2 + qjzero[n];
  569.     }
  570.     return(p/q);
  571. }
  572.  
  573. double pzero(x)
  574. double x;
  575. {
  576. double p, q, z, z2;
  577. int n;
  578.  
  579.     z = 8.0 / x;
  580.     z2 = z * z;
  581.     p = ppzero[5];
  582.     q = qpzero[5];
  583.     for (n=4; n>=0; n--) {
  584.         p = p*z2 + ppzero[n];
  585.         q = q*z2 + qpzero[n];
  586.     }
  587.     return(p/q);
  588. }
  589.  
  590. double qzero(x)
  591. double x;
  592. {
  593. double p, q, z, z2;
  594. int n;
  595.  
  596.     z = 8.0 / x;
  597.     z2 = z * z;
  598.     p = pqzero[5];
  599.     q = qqzero[5];
  600.     for (n=4; n>=0; n--) {
  601.         p = p*z2 + pqzero[n];
  602.         q = q*z2 + qqzero[n];
  603.     }
  604.     return(p/q);
  605. }
  606.  
  607. double yzero(x)
  608. double x;
  609. {
  610. double p, q, x2;
  611. int n;
  612.  
  613.     x2 = x * x;
  614.     p = pyzero[8];
  615.     q = qyzero[8];
  616.     for (n=7; n>=0; n--) {
  617.         p = p*x2 + pyzero[n];
  618.         q = q*x2 + qyzero[n];
  619.     }
  620.     return(p/q);
  621. }
  622.  
  623. double rj0(x)
  624. double x;
  625. {
  626.     if ( x <= 0.0 )
  627.         x = -x;
  628.     if ( x < 8.0 )
  629.         return(jzero(x));
  630.     else
  631.         return( sqrt(TWO_ON_PI/x) *
  632.             (pzero(x)*cos(x-PI_ON_FOUR) - 8.0/x*qzero(x)*sin(x-PI_ON_FOUR)) );
  633.  
  634. }
  635.  
  636. double ry0(x)
  637. double x;
  638. {
  639.     if ( x < 0.0 )
  640.         return(dzero/dzero); /* error */
  641.     if ( x < 8.0 )
  642.         return( yzero(x) + TWO_ON_PI*rj0(x)*log(x) );
  643.     else
  644.         return( sqrt(TWO_ON_PI/x) *
  645.             (pzero(x)*sin(x-PI_ON_FOUR) + 
  646.             (8.0/x)*qzero(x)*cos(x-PI_ON_FOUR)) );
  647.  
  648. }
  649.  
  650.  
  651. double jone(x)
  652. double x;
  653. {
  654. double p, q, x2;
  655. int n;
  656.  
  657.     x2 = x * x;
  658.     p = pjone[8];
  659.     q = qjone[8];
  660.     for (n=7; n>=0; n--) {
  661.         p = p*x2 + pjone[n];
  662.         q = q*x2 + qjone[n];
  663.     }
  664.     return(p/q);
  665. }
  666.  
  667. double pone(x)
  668. double x;
  669. {
  670. double p, q, z, z2;
  671. int n;
  672.  
  673.     z = 8.0 / x;
  674.     z2 = z * z;
  675.     p = ppone[5];
  676.     q = qpone[5];
  677.     for (n=4; n>=0; n--) {
  678.         p = p*z2 + ppone[n];
  679.         q = q*z2 + qpone[n];
  680.     }
  681.     return(p/q);
  682. }
  683.  
  684. double qone(x)
  685. double x;
  686. {
  687. double p, q, z, z2;
  688. int n;
  689.  
  690.     z = 8.0 / x;
  691.     z2 = z * z;
  692.     p = pqone[5];
  693.     q = qqone[5];
  694.     for (n=4; n>=0; n--) {
  695.         p = p*z2 + pqone[n];
  696.         q = q*z2 + qqone[n];
  697.     }
  698.     return(p/q);
  699. }
  700.  
  701. double yone(x)
  702. double x;
  703. {
  704. double p, q, x2;
  705. int n;
  706.  
  707.     x2 = x * x;
  708.     p = 0.0;
  709.     q = qyone[8];
  710.     for (n=7; n>=0; n--) {
  711.         p = p*x2 + pyone[n];
  712.         q = q*x2 + qyone[n];
  713.     }
  714.     return(p/q);
  715. }
  716.  
  717. double rj1(x)
  718. double x;
  719. {
  720. double v,w;
  721.     v = x;
  722.     if ( x < 0.0 )
  723.         x = -x;
  724.     if ( x < 8.0 )
  725.         return(v*jone(x));
  726.     else {
  727.         w = sqrt(TWO_ON_PI/x) *
  728.             (pone(x)*cos(x-THREE_PI_ON_FOUR) - 
  729.                8.0/x*qone(x)*sin(x-THREE_PI_ON_FOUR)) ;
  730.         if (v < 0.0)
  731.             w = -w;
  732.         return( w );
  733.     }
  734. }
  735.  
  736. double ry1(x)
  737. double x;
  738. {
  739.     if ( x <= 0.0 )
  740.         return(dzero/dzero); /* error */
  741.     if ( x < 8.0 )
  742.         return( x*yone(x) + TWO_ON_PI*(rj1(x)*log(x) - 1.0/x) );
  743.     else
  744.         return( sqrt(TWO_ON_PI/x) *
  745.             (pone(x)*sin(x-THREE_PI_ON_FOUR) + 
  746.             (8.0/x)*qone(x)*cos(x-THREE_PI_ON_FOUR)) );
  747. }
  748.  
  749.  
  750. f_besj0()    
  751. {
  752. struct value a;
  753. double x;
  754.     (void) pop(&a);
  755.     if (imag(&a) > zero)
  756.         int_error("can only do bessel functions of reals",NO_CARET);
  757.     push( complex(&a,rj0(real(&a)),0.0) );
  758. }
  759.  
  760.  
  761. f_besj1()    
  762. {
  763. struct value a;
  764. double x;
  765.     (void) pop(&a);
  766.     if (imag(&a) > zero)
  767.         int_error("can only do bessel functions of reals",NO_CARET);
  768.     push( complex(&a,rj1(real(&a)),0.0) );
  769. }
  770.  
  771.  
  772. f_besy0()    
  773. {
  774. struct value a;
  775. double x;
  776.     (void) pop(&a);
  777.     if (imag(&a) > zero)
  778.         int_error("can only do bessel functions of reals",NO_CARET);
  779.     if (real(&a) > 0.0)
  780.         push( complex(&a,ry0(real(&a)),0.0) );
  781.     else {
  782.         push( complex(&a,0.0,0.0) );
  783.         undefined = TRUE ;
  784.     }
  785. }
  786.  
  787.  
  788. f_besy1()    
  789. {
  790. struct value a;
  791. double x;
  792.     (void) pop(&a);
  793.     if (imag(&a) > zero)
  794.         int_error("can only do bessel functions of reals",NO_CARET);
  795.     if (real(&a) > 0.0)
  796.         push( complex(&a,ry1(real(&a)),0.0) );
  797.     else {
  798.         push( complex(&a,0.0,0.0) );
  799.         undefined = TRUE ;
  800.     }
  801. }
  802.  
  803.